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Abstract 

The discrete class algorithm presented in this paper 
is an efficient simulation tool for stochastic processes 
governed by a reasonably small set of transition rates. 
The algorithm is presented, its performance compared 
to prevailing methods and applications to epitaxial 
growth and neuronal models are sketched. 

I. Introduction 

Stochastic processes play a crucial role in many fields 
of science and technology and have received much 
attention ever since the ground-breaking work by 
Einstein, Smoluchowski and others at the beginning 
of the century [1, 2]. While many low-dimensional 
stochastic processes can be treated analytically [3, 4], 
this is no longer the case for spatially extended, 
high-dimensional systems, such as diffusion-limited 
reaction-diffusion systems, epitaxial growth, popula- 
tion dynamics or neuronal interactions [5, 6, 7, 8]. 

The unifying feature of all these systems is that their 
development in time is given by a master equation 

|^(n>*|n (0) ) = 

E W^ n V(n>, t 1 nW) - W a ^ a ,V(n, t | n<°>) 

for the probability of the system to be in state n 
at time t if it was in state n' ' at time t<-°\ Here, 
Wn— >n' are transition rates and n is a vector in an 
m-dimensional discrete state space. Since the mas- 
ter equation can rarely be solved analytically efficient 
simulation methods for the generation of trajectories 
obeying the equation are of tantamount importance. 

In the next section, we present the highly efficient 
discrete class algorithm (DCA) for the simulation of 
systems governed by a reasonably small set of different 



transistion rates. In section 3 we apply the algorithm 
to a simple model of epitaxial growth and demonstrate 
the speed-up compared to prevailing methods. Finally, 
in section 4, we show how the DCA can be used to 
study large neural networks. 

II. The Discrete Class Algorithm 

The discrete class algorithm is an extension of the min- 
imal process method [9] introduced by Gillespie [10]. 
This elegant algorithm proceeds from state n at t to 
n' at t' = t + t as follows: 

1. Calculate the total rate for leaving state n: 

2. Determine the (exponentially distributed) time 
step r = -lnrnd(0, 1]/W n . 

3. Choose a new state n' with probability 
p n (n') = W a ^ n ,/W n . 

The efficiency of this algorithm hinges on the ef- 
ficient implementation of step 3, i.e. the selection of 
the new state n' from the probability distribution 
p n (n') depending on the current state n. For an Tri- 
dimensional problem, the number of states n' acces- 
sible from n (W n ^ n , > 0) will be M ~ O(m). To 
see this, consider a reaction-diffusion system with two 
particle species A and B and one chemical reaction 
A + B — > 0, modelled on a grid of L x L cells. Then, 
the state space is m = 2i 2 -dimensional and as long 
as all cells contain at least one A and one B particle 
each, at every time step one out of M = (4 + 4 + 1) • I? 
possible events has to be chosen: From any one of the 
L 2 cells, either an A or a B particle diffuses to any one 
of the four nearest neighbors or a reaction occurs in it. 
Therefore, linear selection schemes requiring an effort 
of 0(.M) additions per time step are utterly unsuitable 
for large systems, as are rejection methods, which are 



efficient only if p n (n') is restricted to a small inter- 
val [11]. Methods employing binary trees for step 3 
still require an effort of 0(log 2 M) per time step [12]. 

A sophisticated method using a logarithmic classi- 
fication scheme for step 3 has been introduced a few 
years ago [13, 14]. This method yields a computational 
effort for step 3 independent of M, i.e. independent of 
the size of the system, provided the transition rates 
W n _ fn , are independent of M (for a more detailed dis- 
cussion see [11]). The algorithm involves some over- 
head though, and for systems spanning a very large 
range of transition rates (> O(10 14 )), the algorithm 
slows down slightly due to computational precautions 
required to avoid round-off errors. 

The discrete class algorithm is similar in spirit to the 
logarithmic classes, but specifically aimed at systems 
with a reasonably small, discrete set of transition rates, 
i.e. 
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Figure 1: CPU time required per time step for 
different simulation algorithms. All simulations 
were performed on an SGI Indy workstation with 
192 megabyte of memory. 



W n ^ n , e{n,...,r K } 7 K<50 . 

While this may seem a to be a strong restriction at 
first, models of epitaxial growth [15] fulfill these re- 
strictions as well as some neuronal models. 

The DCA implements step 3 of the minimal pro- 
cess method as follows. Each possible transition event 
n — > n' is assigned to one of K classes according to its 
rate 

D v = {n - n' | W n ^ n , = rv} , uGl,...,K . 
Thus the rate of events in class D v is given by 

Ru= E Wn^n' = HA,||r„, 

where \\D V \\ is the number of events in D„, while the 
total rate of events is given by 

K 



Furthermore, within a class each event occurs with 
equal probability. The selection step 3 is thus split 
into two substeps: 

3a. Choose a class D u with probability R v /W a by 
linear selection, i.e. for a p = rnd[0, W n ) select 
that class v for which 



O(K) additions, independent of the size of the system, 
while step 3b requires drawing another uniformly dis- 
tributed random number. Thus, the efficiency of the 
selection algorithm does not depend on the size of the 
system under study. As the number of classes is as- 
sumed to be small, the total rate W n can be calculated 
at every step, keeping numeric inaccuracies to a mini- 
mum. 

III. Comparison with other methods 

To demonstrate the performance of the DCA com- 
pared to the logarithmic classes [5, 11] and other state- 
of-the art methods such as binary trees [12], let us con- 
sider a simple model of epitaxial growth based on [6] : 

• the substrate is an L x L lattice; 

• each lattice site is occupied by one adatom or 
none; 

• in an initial phase, N < L x L adatoms are de- 
posited, which cannot evaporate; 

• an adatom with all four next neighbor sites oc- 
cupied cannot move; 

• all other adatoms diffuse to next neighbor sites 
with rates 
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3b. Select the new state n' from class D v at random. 

The linear selection in step 3a requires drawing 
a single, uniformly distributed random number and 



Here, n is the number of occupied next-neighbor sites, 
h, ks are Planck's and Boltzmann's constants, T is 
temperature and Es, Em are material-dependent en- 
ergies characterizing adatom-substrate and adatom- 
adatom interactions, respectively; both are on the or- 
der of lcV. Thus, after the deposition phase, the sys- 
tem is determined by just five different rates which 



• All input fs , fp consists of delta-spikes, i.e. an 
input event at time T corresponds to the transi- 
t]onv j (T-)->v j {T+) = v j (T-) + l. 

• fp\t) is Poissonian noise with rate l/r p . 

• fs^ (t) is synaptic input from other neurons. 
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Figure 2: Propagation of waves on the retina. 
Grey level indicates time since last firing, black 
being most recent. 



are at T = 600K: w = 3.0 • 10 2 , wx = 1.2 • 1CT 6 , 
w 2 = 4.8 • 1CT 15 , w 3 = 1.9 • 1CT 23 and w 4 = 0. 

Figure 1 shows the CPU time required per step for 
the simulation of this model for different lattice sizes 
using the discrete class, the logarithmic class and the 
binary tree algorithm. This demonstrates clearly the 
superiority of the discrete classes to the other algo- 
rithms in terms of absolute times as well as the size- 
independence of efficiency The minuscule increase in 
CPU time for the DCA at very large lattices is due to 
cache effects, i.e. shortcomings of the hardware; for a 
detailed discussion, see [16]. 

IV. A Neuronal Model 

A crucial problem in modeling the signal processing by 
neuronal networks is the enormous number of neurons 
involved even in simple tasks. Typically, though, only 
a small number of neurons will respond to any one 
stimulus presented e.g. to the eye or the ear. The 
DCA is well suited for the simulation of such largely 
"dormant" systems, since it automatically "focuses" 
on active regions of the system under study. 

To demonstrate the applicability of the DCA to neu- 
ronal studies, we have formulated the following model, 
which is essentially a simplified type of Stein's model 
neuron [17]. 

• At t = 0, each neuron j has the resting mem- 
brane potential Vj(0) = 0. 

• The membrane potential Vj is governed by the 
equation dvj/dt = + fp\t). 



• As the potential reaches a threshold, Vj(t) = 6, 
neuron j fires a spike after an average waiting 
time Tf, which is transmitted to all kj neurons 
receiving input from j; neuron j is reset to an 
absolute refractory state. 

• All input is ignored in the refractory state and 
the neuron returns to the resting state Vj (t) = 
with rate l/r r . 

This model is obviously very well suited for the DCA, 
since it is governed by only three different rates: l/r p , 
1/t/ and l/r r . In order to model spontaneous retinal 
waves as have been observed in newborn ferrets [18], 
we have simulated a grid of neurons with strongly 
localized synaptic connections. The network studied 
had 512 x 512 neurons and some 6.7 • 10 6 synapses, 
the threshold was set to 6 = 7. The simulation was 
stopped after 1.5 million spikes had been generated, 
which required only 80 seconds of CPU time on an 
SGI Indy workstation. Figure 2 shows a typical state 
of activity. 

V. Conclusions 

The DCA algorithm presented here is a powerful tool 
for the study of a large class of stochastic systems and 
should foster research in these fields. The extension of 
the epitaxial model towards more complex phenomena 
is straightforward. 

The neuron model presented above is most likely 
too simplistic to further our understanding of real 
neuronal systems, but we are presently working on 
a faithful implementation of Stein's model. Prelimi- 
nary results indicate that leak currents and inhibitory 
inputs can be included. Inclusion of arbitrary synap- 
tic weights, though, might necessitate recourse to the 
more generally applicable logarithmic class algorithm. 

Note that the effective implementation of the DCA 
requires sophisticated data structures, similar to those 
described in [11]. Source code that can be integrated 
in simulation software via an easy to use interface is 
available from the authors upon request. 
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